<<<<<<< HEAD ======= >>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Wet vs. Dry ======= .bi-alphabet-uppercase::before { content: "\f2a5"; } .bi-alphabet::before { content: "\f68a"; } .bi-amazon::before { content: "\f68d"; } .bi-arrows-collapse-vertical::before { content: "\f690"; } .bi-arrows-expand-vertical::before { content: "\f695"; } .bi-arrows-vertical::before { content: "\f698"; } .bi-arrows::before { content: "\f6a2"; } .bi-ban-fill::before { content: "\f6a3"; } .bi-ban::before { content: "\f6b6"; } .bi-bing::before { content: "\f6c2"; } .bi-cake::before { content: "\f6e0"; } .bi-cake2::before { content: "\f6ed"; } .bi-cookie::before { content: "\f6ee"; } .bi-copy::before { content: "\f759"; } .bi-crosshair::before { content: "\f769"; } .bi-crosshair2::before { content: "\f794"; } .bi-emoji-astonished-fill::before { content: "\f795"; } .bi-emoji-astonished::before { content: "\f79a"; } .bi-emoji-grimace-fill::before { content: "\f79b"; } .bi-emoji-grimace::before { content: "\f7a0"; } .bi-emoji-grin-fill::before { content: "\f7a1"; } .bi-emoji-grin::before { content: "\f7a6"; } .bi-emoji-surprise-fill::before { content: "\f7a7"; } .bi-emoji-surprise::before { content: "\f7ac"; } .bi-emoji-tear-fill::before { content: "\f7ad"; } .bi-emoji-tear::before { content: "\f7b2"; } .bi-envelope-arrow-down-fill::before { content: "\f7b3"; } .bi-envelope-arrow-down::before { content: "\f7b8"; } .bi-envelope-arrow-up-fill::before { content: "\f7b9"; } .bi-envelope-arrow-up::before { content: "\f7be"; } .bi-feather::before { content: "\f7bf"; } .bi-feather2::before { content: "\f7c4"; } .bi-floppy-fill::before { content: "\f7c5"; } .bi-floppy::before { content: "\f7d8"; } .bi-floppy2-fill::before { content: "\f7d9"; } .bi-floppy2::before { content: "\f7e4"; } .bi-gitlab::before { content: "\f7e5"; } .bi-highlighter::before { content: "\f7f8"; } .bi-marker-tip::before { content: "\f802"; } .bi-nvme-fill::before { content: "\f803"; } .bi-nvme::before { content: "\f80c"; } .bi-opencollective::before { content: "\f80d"; } .bi-pci-card-network::before { content: "\f8cd"; } .bi-pci-card-sound::before { content: "\f8ce"; } .bi-radar::before { content: "\f8cf"; } .bi-send-arrow-down-fill::before { content: "\f8d0"; } .bi-send-arrow-down::before { content: "\f8d1"; } .bi-send-arrow-up-fill::before { content: "\f8d2"; } .bi-send-arrow-up::before { content: "\f8d3"; } .bi-sim-slash-fill::before { content: "\f8d4"; } .bi-sim-slash::before { content: "\f8d5"; } .bi-sourceforge::before { content: "\f8d6"; } .bi-substack::before { content: "\f8d7"; } .bi-threads-fill::before { content: "\f8d8"; } .bi-threads::before { content: "\f8d9"; } .bi-transparency::before { content: "\f8da"; } .bi-twitter-x::before { content: "\f8db"; } .bi-type-h4::before { content: "\f8dc"; } .bi-type-h5::before { content: "\f8dd"; } .bi-type-h6::before { content: "\f8de"; } .bi-backpack-fill::before { content: "\f8df"; } .bi-backpack::before { content: "\f8e0"; } .bi-backpack2-fill::before { content: "\f8e1"; } .bi-backpack2::before { content: "\f8e2"; } .bi-backpack3-fill::before { content: "\f8e3"; } .bi-backpack3::before { content: "\f8e4"; } .bi-backpack4-fill::before { content: "\f8e5"; } .bi-backpack4::before { content: "\f8e6"; } .bi-brilliance::before { content: "\f8e7"; } .bi-cake-fill::before { content: "\f8e8"; } .bi-cake2-fill::before { content: "\f8e9"; } .bi-duffle-fill::before { content: "\f8ea"; } .bi-duffle::before { content: "\f8eb"; } .bi-exposure::before { content: "\f8ec"; } .bi-gender-neuter::before { content: "\f8ed"; } .bi-highlights::before { content: "\f8ee"; } .bi-luggage-fill::before { content: "\f8ef"; } .bi-luggage::before { content: "\f8f0"; } .bi-mailbox-flag::before { content: "\f8f1"; } .bi-mailbox2-flag::before { content: "\f8f2"; } .bi-noise-reduction::before { content: "\f8f3"; } .bi-passport-fill::before { content: "\f8f4"; } .bi-passport::before { content: "\f8f5"; } .bi-person-arms-up::before { content: "\f8f6"; } .bi-person-raised-hand::before { content: "\f8f7"; } .bi-person-standing-dress::before { content: "\f8f8"; } .bi-person-standing::before { content: "\f8f9"; } .bi-person-walking::before { content: "\f8fa"; } .bi-person-wheelchair::before { content: "\f8fb"; } .bi-shadows::before { content: "\f8fc"; } .bi-suitcase-fill::before { content: "\f8fd"; } .bi-suitcase-lg-fill::before { content: "\f8fe"; } .bi-suitcase-lg::before { content: "\f8ff"; } .bi-suitcase::before { content: "\f900"; } .bi-suitcase2-fill::before { content: "\f901"; } .bi-suitcase2::before { content: "\f902"; } .bi-vignette::before { content: "\f903"; } >>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843

Wet vs. Dry

<<<<<<< HEAD
======= >>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843

Comparing various definitions of a “wet sample”.

Load things

Packages, spatial layers, etc.

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
library(tidyverse)
library(sf)
library(mapview)
library(tbeptools)
library(patchwork)

dbasins <- read_sf(here::here("data", "geospatial", "dbasins-galb.shp"))
catchments <- read_sf(here::here("data", "geospatial", "TB_Select_Fib_Catchments_Dissolved.shp"))
swfwmd_pixels <- read_sf(here::here("data", "geospatial", "swfwmd_pixel_2_utm_m_83.shp"))

Source functions

These will eventually be incorporated into tbeptools; for now sourcing the .R files

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# to download Entero data at specified stations
source(here::here("R", "read_importwqpStationsAndStuff.R"))

# to download precip data from SWFWMD
source(here::here("R", "raindl.R"))

# to identify wet/dry fib samples  
source(here::here("R", "anlz_fibwetdry.R"))

Download things

Enterococcus data

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
stns <- unique(catch_pixels$station)
stns[which(stns == "21FLCOSP_WQX-CENTRAL *")] <- "21FLCOSP_WQX-CENTRAL CANAL"
# define arguments


stations <- stns  # can subset if desired
entero_names <- c("Enterococci",
                  "Enterococcus")
startDate <- as.Date("2018-01-01")
endDate <- as.Date("2022-12-31")
<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# put them together
# date format has to be mm-dd-yyyy for the url
args <- list(
    siteid = stations,
    characteristicName = entero_names,
    startDateLo = format(startDate, "%m-%d-%Y"),
    startDateHi = format(endDate, "%m-%d-%Y")
)
<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# grab the data
fibdat <- read_importwqpStationsAndStuff(args = args)

# fibdat2 <- read_importwqpStationsAndStuff(args = args)

fibdat$date <- lubridate::ymd(fibdat$date)
<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
fibdat2 <- fibdat |> 
    mutate(year = lubridate::year(date))
table(fibdat2$station, fibdat2$year)

Precip data

Actually need to loop through some years

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# yrs <- 2018:2022
# 
# prcp_out <- list()
# 
# for(i in seq_along(yrs)){
#     yr = yrs[i]
#     prcptmp <- read_importrain(curyr = yr, catch_pixels = catch_pixels,
#                            quiet = FALSE)
#     prcp_out[[i]] <- prcptmp
# }
# 
# prcpdat <- bind_rows(prcp_out)
# 
# save(prcpdat, file = here::here("data", "precipitation", "daily_rain_stns_2018-2022.RData"))
load(here::here("data", "precipitation", "daily_rain_stns_2018-2022.RData"))

Fix Central Canal in the precip data station listing, so it will match the station in the fib data.

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
prcpdat <- prcpdat |> 
    mutate(station = case_when(station == "21FLCOSP_WQX-CENTRAL *" ~ "21FLCOSP_WQX-CENTRAL CANAL",
                               .default = station))

Define stations

Some of the FIB stations are in the same catchments used for rain data - when calculating FIB amounts, we want them included; but when calculating rain amounts, we don’t want to falsely replicate them. So define which ones should be removed for calculating rain stats:

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
stn_prcp_dupes <- c(
    '21FLHILL_WQX-103',
    '21FLHILL_WQX-619',
    '21FLHILL_WQX-597',
    '21FLHILL_WQX-264', 
    '21FLHILL_WQX-180', 
    '21FLDOH_WQX-MANATEE153',
    '21FLCOSP_WQX-SOUTH CANAL',
    '21FLCOSP_WQX-CENTRAL CANAL',
    '21FLPDEM_WQX-39-05',
    '21FLMANA_WQX-LM4', 
    '21FLMANA_WQX-LM5',
    '21FLMANA_WQX-LM6'
)

Create combinations of definitions

Remember, temporal threshold of 1 means “sample day” - for sample day + previous day, use 2 for the threshold. Previous 2 days, 3. Previous 3 days, 4.

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
definition_combos <- expand.grid(wet_threshold = c(0.25, 0.5, 0.75, 1), 
                                 temporal_threshold = 1:4)

Loop through them

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
outs_df <- list()
outs_comp_df <- list()
outs_graphics_conc <- list()
outs_graphics_prcp <- list()
for(i in 1:nrow(definition_combos)){
    wetthrs <- definition_combos$wet_threshold[i]
    tmpthrs <- definition_combos$temporal_threshold[i]
    
    fib_out <- anlz_fibwetdry(fibdat, prcpdat, 
                        temporal_window = tmpthrs, 
                        wet_threshold = wetthrs)
    
    
    # sample information
    # ------------------------------------------------------------------------
    # number wet samples
    nwet <- sum(fib_out$wet_sample, na.rm = TRUE)
    # number dry samples
    ndry <- sum(fib_out$wet_sample == FALSE, na.rm = TRUE)
    # number samples
    nsamps <- nrow(fib_out)
    # proportion of samples that are wet, out of those that could be ID'd:
    prop_wet <- nwet / (nwet + ndry)
    # percent wet
    pct_wet <- round(prop_wet*100, 1)
    # number stations
    nstns <- length(unique(fib_out$station))
    
    overall <- data.frame(
        wet_threshold = wetthrs,
        temporal_window = tmpthrs,
        nsamps,
        nstns,
        nwet,
        ndry,
        prop_wet,
        pct_wet
    )

    
    # concentration info
    # -------------------------------------------------------------------------
    conc_out <- fib_out |> 
        filter(!is.na(wet_sample),
               !is.na(FIBconc)) |> 
        summarize(.by = wet_sample,
                  nsamples = n(),
                  medianEntero = median(FIBconc),
                  IQREntero = IQR(FIBconc),
                  pct75Entero = quantile(FIBconc, 0.75),
                  geom_meanEntero = exp(mean(log(FIBconc + 0.001))),
                  quantilesEntero = list(quantile(FIBconc)),
                  nexceedances = sum(FIBconc > 70),
                  propExceedances = nexceedances / nsamples)
    
    
    # precip info
    # --------------------------------------------------------------------------
    
    # remove dupe stations, so don't double-count rain from same catchments
    prcp_out <- fib_out |> 
        filter(!(station %in% stn_prcp_dupes),
               !is.na(wet_sample)) |>
        select(rain_sampleDay, rain_total, wet_sample) |> 
        summarize(.by = wet_sample,
                  median_prcp = median(rain_total),
                  quantiles_prcp = list(quantile(rain_total)))
    
    
    # combine for writing out  
    comb_out <- full_join(conc_out, prcp_out,
                          by = "wet_sample") |> 
        mutate(wet_threshold = wetthrs,
               temporal_window = tmpthrs) |> 
        relocate(c(wet_threshold, temporal_window))
    
    outs_comp_df[[i]] <- comb_out
    
    
    comp2 <- conc_out |> 
        select(-IQREntero, -quantilesEntero, -nexceedances) |> 
        pivot_longer(-wet_sample,
                     names_to = "stat",
                     values_to = "value") |>
        mutate(wet_sample = case_when(wet_sample == TRUE ~ "wet",
                                      wet_sample == FALSE ~ "dry",
                                      .default = "somethingwrong")) |> 
        pivot_wider(names_from = c(stat, wet_sample),
                    values_from = value)
    
    outs_df[[i]] <- cbind(overall, comp2)

    
    # graphs
    # --------------------------------------------------------------------------
    
    outs_graphics_conc[[i]] <- ggplot(fib_out, aes(x = wet_sample,
                                                   y = FIBconc,
                                                   fill = wet_sample)) +
        geom_violin(draw_quantiles = c(0.25, 0.5, 0.75),
                    alpha = 0.3) +
        scale_y_log10() +
        labs(subtitle = glue::glue("\n{pct_wet}% of {nsamps} samples across {nstns} stations 'wet' \n'wet' threshold {wetthrs} inches over {tmpthrs} days")) +
        theme(legend.position = "none")
    
    
    outs_graphics_prcp[[i]] <- ggplot(fib_out, aes(x = wet_sample,
                                                   y = rain_total)) +
        geom_violin(draw_quantiles = c(0.25, 0.5, 0.75)) +
        labs(subtitle = glue::glue("'wet' threshold {wetthrs} inches over {tmpthrs} days"))
}
<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
wet_drys <- bind_rows(outs_comp_df) |> 
    arrange(wet_threshold,
            temporal_window,
            wet_sample) |> 
    left_join(bind_rows(outs_df))


overall <- bind_rows(outs_df) |> 
    mutate(threshold = glue::glue("{wet_threshold} inches over {temporal_window} days"),
           median_diff = medianEntero_wet - medianEntero_dry,
           geomMean_diff = geom_meanEntero_wet - geom_meanEntero_dry,
           pct75_diff = pct75Entero_wet - pct75Entero_dry,
           excd_diff = propExceedances_wet - propExceedances_dry)
<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
overall |> 
    select(threshold,
           ends_with("diff"),
           everything()) |> 
    View()

Graphs!

Visualize Entero concentrations under various thresholds

These are violin plots - like boxplots, but showing the shape of the distribution - with quantiles drawn as horizontal lines (25th %ile; median; 75th %ile)

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
grps <- outs_graphics_conc

(grps[[1]] + grps[[2]]) / (grps[[3]] + grps[[4]])  +
    plot_annotation(title = "varying precip over 1 day")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
(grps[[5]] + grps[[6]]) / (grps[[7]] + grps[[8]]) +
    plot_annotation(title = "varying precip over 2 days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
(grps[[9]] + grps[[10]]) / (grps[[11]] + grps[[12]]) +
    plot_annotation(title = "varying precip over 3 days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
(grps[[13]] + grps[[14]]) / (grps[[15]] + grps[[16]]) +
    plot_annotation(title = "varying precip over 4 days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843

It looks like the shape doesn’t change too much by varying the amount of precip over the same number of days. Let’s look at a half-inch threshold over each of the days:

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
(grps[[2]] + grps[[6]]) / (grps[[10]] + grps[[14]]) +
    plot_annotation(title = "0.5 inches over varying # of days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843

It’s hard to compare these because of the NAs for wet-sample throwing off the x-axis. I’ve kept those because it’s important to notice that we’re getting NAs when we do this - presumably only at the beginning of the entire time series when we don’t have data from the days before sampling because of the year cutoff. Important to note that we need to download precip data for at least several days before FIB samples - e.g. if you want to look at FIBs from 2018-2023, download precip data starting in mid-December 2017, to avoid those NAs.

Removing unidentified samples

Make those graphs again but cutting out the NAs now.

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
outs_graphics_conc2 <- list()
for(i in 1:nrow(definition_combos)){
    wetthrs <- definition_combos$wet_threshold[i]
    tmpthrs <- definition_combos$temporal_threshold[i]
    
    fib_out <- anlz_fibwetdry(fibdat, prcpdat, 
                        temporal_window = tmpthrs, 
                        wet_threshold = wetthrs) |> 
        filter(!is.na(wet_sample))
    
    
    # sample information
    # ------------------------------------------------------------------------
    # number wet samples
    nwet <- sum(fib_out$wet_sample, na.rm = TRUE)
    # number dry samples
    ndry <- sum(fib_out$wet_sample == FALSE, na.rm = TRUE)
    # number samples
    nsamps <- nrow(fib_out)
    # proportion of samples that are wet, out of those that could be ID'd:
    prop_wet <- nwet / (nwet + ndry)
    # percent wet
    pct_wet <- round(prop_wet*100, 1)
    # number stations
    nstns <- length(unique(fib_out$station))
    
    
    # graphs
    # --------------------------------------------------------------------------
    
    outs_graphics_conc2[[i]] <- ggplot(fib_out, aes(x = wet_sample,
                                                   y = FIBconc,
                                                   fill = wet_sample)) +
        geom_violin(draw_quantiles = c(0.25, 0.5, 0.75),
                    alpha = 0.3) +
        scale_y_log10() +
        labs(subtitle = glue::glue("\n{pct_wet}% of {nsamps} samples across {nstns} stations 'wet' \n'wet' threshold {wetthrs} inches over {tmpthrs} days")) +
        theme(legend.position = "none")

}
<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
grps <- outs_graphics_conc2

# over one day
(grps[[1]] + grps[[2]]) / (grps[[3]] + grps[[4]]) +
    plot_annotation(title = "varying precip over 1 day")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# 2 days
(grps[[5]] + grps[[6]]) / (grps[[7]] + grps[[8]]) +
    plot_annotation(title = "varying precip over 2 days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# 3 days
(grps[[9]] + grps[[10]]) / (grps[[11]] + grps[[12]]) +
    plot_annotation(title = "varying precip over 3 days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# 4 days    
(grps[[13]] + grps[[14]]) / (grps[[15]] + grps[[16]]) +
    plot_annotation(title = "varying precip over 4 days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# 0.5 in
(grps[[2]] + grps[[6]]) / (grps[[10]] + grps[[14]]) +
    plot_annotation(title = "0.5 inches over varying # of days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
# 0.75 in
(grps[[3]] + grps[[7]]) / (grps[[11]] + grps[[15]]) +
    plot_annotation(title = "0.75 inches over varying # of days")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843

Distribution of concentrations really doesn’t change too much, visually. It’s the top end of wet samples that changes shape a bit, and the percent of samples that get classified as wet.

~18% of samples get classified as ‘wet’ with 0.5 inch over 2 days or 0.75 inches over 3 days.

~25% of samples are classified as ‘wet’ with 0.5 inch over 3 days.

These numbers all feel pretty reasonable.

Proportion of Exceedances?

Does this change?

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
long <- overall |> 
    select(threshold,
           wet_threshold,
           temporal_window,
           propExceedances_dry,
           propExceedances_wet,
           geommean_dry = geom_meanEntero_dry,
           geommean_wet = geom_meanEntero_wet) |> 
    pivot_longer(-c(threshold, wet_threshold, temporal_window),
                 names_to = c("param", "sampleType"),
                 names_sep = "_",
                 values_to = "value") |> 
    pivot_wider(names_from = param,
                values_from = value)
<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
ggplot(long, aes(x = sampleType,
                 y = propExceedances,
                 col = as.factor(wet_threshold),
                 shape = as.factor(temporal_window))) +
    geom_point(size = 3,
               alpha = 0.8) +
    labs(title = "Exceedances",
         x = "Sample Type",
         y = "proportion of samples over 70 CFU/100mL",
         shape = "# days",
         col = "inches of precip")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
ggplot(filter(long,
              wet_threshold %in% c(0.5, 0.75),
              temporal_window %in% c(2, 3)), 
       aes(x = sampleType,
                 y = propExceedances,
                 col = as.factor(wet_threshold),
                 shape = as.factor(temporal_window))) +
    geom_point(size = 3,
               alpha = 0.8) +
    labs(title = "Exceedances",
         subtitle = "only 0.5 or 0.75 inches",
         x = "Sample Type",
         y = "proportion of samples over 70 CFU/100mL",
         shape = "# days",
         col = "inches of precip")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843

The biggest jumps of exceedance seem to happen by increasing the number of inches required to classify a sample as ‘wet’. The difference seems to be around 10% between 4 days and 1 day within each “inches” threshold. Whichever threshold combinations get used, there will be more exceedances in wet samples than dry.

Half an inch over 2 days seems essentially the same across all these metrics as 0.75 inches over 3 days.

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
(grps[[6]] + grps[[11]]) / (grps[[10]] + plot_spacer()) +
    plot_annotation(title = "Comparison of 0.5 / 0.75 inches")
<<<<<<< HEAD

=======

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843

Let’s examine any other differences via table (e.g. geometric mean)

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
overall |> 
    filter(wet_threshold %in% c(0.5, 0.75),
              temporal_window %in% c(2, 3)) |> 
    select(threshold, everything()) |> 
    DT::datatable()
<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843

Again, 0.5 inches over 2 days is pretty similar to 0.75 inches over 3 days. The others have pretty different geometric means and number of samples classified as wet.

We need a certain amount of samples to even be considered wet before we can reliably distinguish - e.g. if only 5% of samples are wet, we can’t really say much. So what have we got? The histogram shows proportions (e.g. 0.1 = 10%).

<<<<<<< HEAD
=======
>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843 Code 
hist(overall$prop_wet, breaks = 10)
<<<<<<< HEAD

=======

Eval a large grid of options

Code 
wetrng <- c(0.25, 4)
wetstp <- 0.25
dys <- c(1, 14)
definition_combos <- expand.grid(wet_threshold = seq(wetrng[1], wetrng[2], by = wetstp), 
                                 temporal_threshold = dys[1]:dys[2])

These examples show a large range of options for number of preceding days and cumulative precipitation to assess the “response surface” for samples defined as wet or dry and their differences. The goal is to better understand how the selected options can discriminate between samples or to determine when the selected options are not useful. Every combination of options from 0.25 to 4 inches of precipitation increasing by 0.25 inches over 1 to 14 days are evaluated for a total of 224 combinations.

Code 
outs_df <- list()
for(i in 1:nrow(definition_combos)){
    # cat(i, '\t')
    wetthrs <- definition_combos$wet_threshold[i]
    tmpthrs <- definition_combos$temporal_threshold[i]
    
    fib_out <- anlz_fibwetdry(fibdat, prcpdat, 
                              temporal_window = tmpthrs, 
                              wet_threshold = wetthrs)
    
    
    # sample information
    # ------------------------------------------------------------------------
    # number wet samples
    nwet <- sum(fib_out$wet_sample, na.rm = TRUE)
    # number dry samples
    ndry <- sum(fib_out$wet_sample == FALSE, na.rm = TRUE)
    # number samples
    nsamps <- nrow(fib_out)
    # proportion of samples that are wet, out of those that could be ID'd:
    prop_wet <- nwet / (nwet + ndry)
    # percent wet
    pct_wet <- prop_wet*100
    # number stations
    nstns <- length(unique(fib_out$station))
    
    overall <- data.frame(
        wet_threshold = wetthrs,
        temporal_window = tmpthrs,
        nsamps,
        nstns,
        nwet,
        ndry,
        prop_wet,
        pct_wet
    )
    
    
    # concentration info
    # -------------------------------------------------------------------------
    conc_out <- fib_out |> 
        filter(!is.na(wet_sample),
               !is.na(FIBconc)) |> 
        summarize(.by = wet_sample,
                  nsamples = n(),
                  medianEntero = median(FIBconc),
                  IQREntero = IQR(FIBconc),
                  pct75Entero = quantile(FIBconc, 0.75),
                  geom_meanEntero = exp(mean(log(FIBconc + 0.001))),
                  quantilesEntero = list(quantile(FIBconc)),
                  nexceedances = sum(FIBconc > 70),
                  propExceedances = nexceedances / nsamples)
    
    
    # precip info
    # -------------------------------------------------------------------------
    
    
    comp2 <- conc_out |> 
        select(-IQREntero, -quantilesEntero, -nexceedances) |> 
        pivot_longer(-wet_sample,
                     names_to = "stat",
                     values_to = "value") |>
        mutate(wet_sample = case_when(wet_sample == TRUE ~ "wet",
                                      wet_sample == FALSE ~ "dry",
                                      .default = "somethingwrong")) |> 
        pivot_wider(names_from = c(stat, wet_sample),
                    values_from = value)
    
    if(sum(names(comp2) %in% c('nsamples_dry', 'nsamples_wet')) == 2)
        comp2 <- comp2 |> 
            mutate(
                meddiff = medianEntero_wet - medianEntero_dry,
            )
        
    outs_df[[i]] <- cbind(overall, comp2)
    
    
}

# combined results
res <- do.call('bind_rows', outs_df) |> 
    mutate(
        pct_wet = ifelse(pct_wet == 0, NA, pct_wet), 
        pct_lab = case_when(
            round(pct_wet, 0) == 0 ~ '< 1%', 
            is.na(pct_wet) ~ '', 
            T ~ paste0(round(pct_wet, 0), '%')
            )
    )
Code 
ggplot(res, aes(x = factor(wet_threshold), y = factor(temporal_window))) +
    geom_tile(aes(fill = pct_wet), color = 'black') +
    geom_text(aes(label = pct_lab), size = 3, color = 'white') +
    scale_fill_viridis_c() +
    scale_x_discrete(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    theme(legend.position  = 'top') + 
    labs(
        y = "Temporal Window (days)",
        x = "Wet Threshold (inches)",
        fill = 'Percent samples "Wet"'
    )

Code 
ggplot(res, aes(x = factor(wet_threshold), y = factor(temporal_window))) +
    geom_tile(aes(fill = meddiff), color = 'black') +
    geom_text(aes(label = round(meddiff, 0)), size = 3, color = 'white') +
    scale_fill_viridis_c() +
    scale_x_discrete(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    theme(legend.position  = 'top') + 
    labs(
        y = "Temporal Window (days)",
        x = "Wet Threshold (inches)",
        fill = 'Diff. in median between wet/dry'
    )

Code 
ggplot(res, aes(x = factor(wet_threshold), y = factor(temporal_window))) +
    geom_tile(aes(fill = propExceedances_wet), color = 'black') +
    geom_text(aes(label = round(propExceedances_wet, 2)), size = 3, color = 'white') +
    scale_fill_viridis_c() +
    scale_x_discrete(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    theme(legend.position  = 'top') + 
    labs(
        y = "Temporal Window (days)",
        x = "Wet Threshold (inches)",
        fill = 'Prop. exceedances in wet samples'
    )

Code 
toplo <- res |> 
    select(wet_threshold, temporal_window, pct_wet, meddiff, propExceedances_wet) |>
    pivot_longer(-c(wet_threshold, temporal_window),
                 names_to = "stat",
                 values_to = "value") |> 
    mutate(
        stat = factor(stat, levels = c('pct_wet', 'meddiff', 'propExceedances_wet'), 
                      labels = c('Percent wet samples', 'Diff. in median between wet/dry', 'Prop. exceedances in wet samples'))
    )

ggplot(toplo, aes(x = wet_threshold, y = value, group = factor(temporal_window), color = temporal_window)) +
    geom_line() +
    scale_color_viridis_c() +
    theme_minimal() +
    theme(
        legend.position  = 'top', 
        strip.placement = 'outside'
        ) + 
    facet_wrap(~stat, scales = 'free_y', strip.position = 'left') +
    labs(
        y = NULL,
        x = "Wet Threshold (inches)",
        color = 'Temporal Window (days)'
    )

These examples show competing objectives for defining wet/dry samples, e.g., maximizing median differences while creating a balance between wet/dry samples. For example, the maximum medium difference has a very small percentage of wet samples because the criteria for wet samples is very strict, e.g., an excessive amount of rain in a very short period of time. The criteria can be optimized to identify which combination of wet threshold and temporal window maintains the best balance in percentage between wet/dry samples and the greatest median difference.

>>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843
<<<<<<< HEAD ======= >>>>>>> 73317b75bfa7d009ad3711bcd7401729aca24843